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We present a fast, arbitrarily accurate method to simulate the effect of gravitational lensing of 
the Cosmic Microwave Background anisotropies and polarization fields by large scale structures. 
We demonstrate the efficiency and accuracy of the method and exhibit their dependence on the 
algorithm parameters. 
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I. INTRODUCTION 

Weak lensing effects on the Cosmic Microwave 
I Background (CMB) temperature and polarization 
anisotropies has been proposed as a probe of the total 
[ matter distribution in Large Scale Structures (LSS) be- 
tween us and the surface of last scattering [l|, 0, H, [1, d, 

; S0JJ,[i3,[il,[il,[i3,[il,[il,[il,[ii,[23,M- Although 

sensitive to the cumulative distribution of matter, it is 
quite complementary to the other probes of the matter 

I distribution of the LSS. Indeed, it does not suffer from 
bias effects (as e.g. galaxy redshift surveys, Lyman-a for- 

I est), or from possible mis-determination of the redshift 
sources (cosmic shear measurements on galaxies) . In ad- 
dition, due to the high redshift of the source (last scatter- 
ing surface) and the lensing efficiency function, weak lens- 
ing of CMB anisotropies is mostly sensitive to large scale 
structures which are still (mainly) in the linear regime, 
which makes it a very useful tool for cosmology, in par- 

' ticular to constrain the properties of neutrinos [2 2 . |23| . 
Unlike shear measurements on galaxies, where the (re- 

i duced) shear field is directly sampled via measurement 
of galaxy ellipticities, measuring weak lensing effects on 

] the CMB is complicated by the fact that the source it- 
self can only be modeled by a stochastic realization of 
a field. However, theoretical arguments lead us to think 
that CMB anisotropies are highly Gaussian [23, |2^ [2^ , 
which has been confirmed on the data, at large scales, 
using different non-Gaussianity estimators (pdf, bispcc- 
trum, wavelet skewness and kurtosis, Minkowski func- 

[ tionals, etc.). These properties of the CMB anisotropies 
can be used to disentangle to some extent the stochastic 
properties of the (unlensed) CMB anisotropies from the 
stochastic properties of the lens (i.e. the LSS) as the lens- 
ing effect induces small specific non-Gaussian features in 
the CMB maps: it locally correlates the anisotropies with 
their gradient 0, [13, [ll, [H, [s^l which in turn has lead 
to the development of specific estimators of the lensing 
potential field and its power spectrum [ll|, [l^, [ij, [3l| . 

Recently, weak lensing of the CMB anisotropics by 
LSS has been measured on WMAP data [2l| by cross- 
correlation with a high-redshift radio galaxy catalog. Al- 
though marginally detectable in WMAP data due to its 
noise level, CMB lensing should be measured with high 
signal to noise by Planck with temperature anisotropies 
[Hi . [3^ without needing to rely on an external data 



set. However, in order to carry out such a measure- 
ment on realistic CMB data, the impact of instrumental 
(anisotropic beams, missing data, correlated noise) and 
astrophysical (Galaxy contamination, point sources, etc.) 
systematic effects on the CMB lensing estimators has to 
be studied with great care. 

The power spectra (temperature, and polarization) can 
be computed, using simple Taylor expansion at large 
scales [lO[, or a more clever resummation scheme at 
smaller scales where the displacement field amplitude is 
comparable to the wavelength of the anisotropies [l^. 
For smaller scales, or to investigate the different system- 
atics described above, the development of fast and ac- 
curate methods to simulate the lensed CMB maps are 
needed. 

This simulation is two-fold. On one hand, an accurate 
simulation of the large scale structure induced lensing de- 
flection field is needed. On the other hand, one needs a 
method to apply this deflection fleld to an unlensed, sim- 
ulated, CMB. We will not consider the first part of this 
program. Indeed, approximating the lensing effect with 
a single lens plane in the so-called Born approximation 
[l0| has been shown to be an excellent approximation, 
both for temperature and polarization anisotropics. In 
this case, the simulation of lensed CMB maps reduces to 
an accurate resampling of the unlensed anisotropies at 
displaced positions. To solve this last problem, several 
technical solutions have been implemented. In the pub- 
licly distributed Lenspix code [33|, different possibilities 
are available, namely: 

• brute-force resampling by direct resummation of 
spherical harmonics at displaced positions (slow, 
but very accurate, this option should be consid- 
ered as the "benchmark" for all other resampling 
methods) 

• resampling on locally Cartesian grids with subse- 
quent polynomial interpolation 

For the last option, an interesting speed-up has been pro- 
posed by Hirata [3, [l3| by noting that a band-limited 
signal in spherical harmonics can be recast as a band- 
limited signal in regular Fourier modes on a {9,ip), thus 
allowing a fast resampling of the signal on a Cartesian 
(6*, if) grid using 2D FFTs. 

In this paper, we investigate a variation on Hirata's 
idea (l6l . \u\ , where the oversampling plus polynomial in- 
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tcrpolation is replaced by an approximate (but arbitrarily 
accurate) Fast Fourier Transform (FFT) resampling on 
irregularly spaced grid points . 

The rest of this paper is organized as follows: In section 
|TI]we briefly describe the resampling technique (hereafter 
NFFT) . This is followed by a brief description of the weak 
lensing of primary CMB fields in section IIIII We also 
describe how remapping of CMB fields on the surface of 
the unit sphere can be recast into remapping of the latter 
on the surface of a 2-d torus. Section IIVI describes the 
details of the simulation procedure for lensed CMB fields 
using NFFT on the surface of a 2-d torus. Finally we 
summarize our results in section |Vl 



II. NON-EQUISPACED FAST FOURIER 
TRANSFORM (NFFT) 

The fast Fourier transform for non-equispaced grid 
points (NFFT) is a generahzation of FFT [11, [Hi. The 
essential idea is that of approximating the reproducing 
kernel of the standard FFT [s^ using a window func- 
tion of specific properties. Suppose we know a function 
/ through N evaluation fk in the frequency domain. Ac- 
cording to NFFT, Fourier transform of that function eval- 
uated at M non-equispaced grid points in spatial domain 
can be written as, 



turns out to be the best. It has been shown that for a 
fixed oversampling factor a > 1, the approximation error 
decays exponentially with K [ssl . [36l | . 

III. WEAK LENSING OF CMB 

The CMB radiation field is completely characterized 
by its temperature anisotropy, T{0, ip) , and polarization, 
P{0,(p), in each direction on the sky. Since temperature 
anisotropy is a spin-0 field on the sphere, it can be con- 
veniently expanded in spin-0 spherical harmonics, 

Tie, if) = J2 TlraYlm{e,f) (2) 

The polarization field can be described by the Stokes 
parameters, Q{d,ip) and U{9,ip), with respect to a par- 
ticular choice of coordinate system on the sky. One 
can conveniently combine the Stokes parameters into 
a single complex quantity representing the polariza- 
tion, P{9,Lp) = (Q + iU) {9 , (f) . Due to its transfor- 
mation properties under rotations, the polarization P 
is a spin-2 field on the sphere. One may thus expand 
P(9. ip) in terms of spin-2 spherical harmonics, 2Yi,n{9-, f) 
[liSliOl, as 
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Here the window function (f>{^) has compact support 
[—a, a] and its Fourier transform 4>{x) assumes small val- 
ues outside some interval [—K, K]. a is the over-sampling 
factor and it is required to avoid the aliasing error. A 
convenient choice for cr is 2, however cr = 3/2 is sufficient 
to get good accuracy, a has to be chosen slightly smaller 
than TT (2 — 2 /a). Since the evaluation of the summation 
over k requires an equispaced FFT of length aN, (j){^) 
has to be well localized in k-space in order to avoid the 
aliasing error with minimal computational cost. On the 
other hand, the summation over m can be evaluated with 
minimum truncation error if the window function is well 
localized in the spatial domain. Hence, the efficient eval- 
uation of f{x) on irregularly spaced grid points requires 
a window function that is well localized in both space 
and frequency domain. It has computational complexity 
O {aN log N + K M) where K is the number of terms 
considered in the spatial approximation, M is the num- 
ber of real space samples, and N the number of Fourier 
modes. Among a number of window functions (Gaus- 
sian, B-spline, Sinc-power, Kaiscr-Bessel), Kaiscr-Bessel 



Pi9,ip) = iQ + iU)i9,ip) 



(3) 
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In the above equation, 2Pim = ^{Eim + iBi„i), where 
Ei„i and Bi„i are the electric and magnetic modes of the 
polarization field in harmonic space. 

Weak lensing induces a deflection field d{9,ip), i.e. a 
mapping between the direction of a given light ray on 
the last scattering surface and the direction in which we 
observe it. Since the deflection field is a vector field on 
the sphere, it can be decomposed in terms of gradient- 
free and curl-free components in the most general form 
as, 



(4) 
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where $(0, (p) and x{9, cp) are two scalar fields on the 
sphere. £a b is the covariant antisymmetric tensor of rank 
2 on the unit sphere. In terms of null basis vectors (m, ffi) 
which define a diad on the unit sphere, eat can be ex- 
pressed as. 



,6 = i{mamb - niamb) 



(5) 



The gradient-free com pon ent can be ignored as it is 
negligible in most cases [30| and is exactly zero in the 
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Born approximation that we use here, as this term can 
only arise when taking into account the lens-lens cou- 
plings. In the Born approximation, the Icnsing deflec- 
tion is calculated on the unlcnsed line of sight so the 
lensed map is a local function of the deflection vector, 
da{d,(p) = Va$(6', v)i where $(0, (p) is the lensing po- 
tential. This projected $(6*, Lp) potential is related to the 

3-d the gravitational potential '^[D, D{D, 6, ip)) as, 
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determined by the deflection field d{6^ (p), 

cos 9' = cos c? cos — sin d sin 6* cos a (11) 
sin a sin((i) 



smf 



(12) 



The extra factor exp [— 2i(7 — a)], that appears in case 
of polarization [lOj, is there to rotate the basis vectors 
{eg', e^') at {9' , cj)') to match them with the basis vectors 
ieg,e^) at {9,(p). 



^iD,D{D,9,ip)) (6) 



where D is the comoving coordinate distance along the 
line of sight and Da is the comoving angular diameter 
distance associated with D. Ds is the coordinate distance 
to the last scattering surface. 

Similarly to CMB temperature anisotropy, the Icnsing 
potential transforms like a spin-zero field on the sphere. 
Hence it may also be expanded in spin-0 spherical har- 
monics. 



Since the deflection field d{6, ip) is a vector field on the 
sphere, it can be expanded in spin-1 spherical harmonics. 



da{9,^) = Va$(^,V5) 



1=0 m=-l 

X [(-i)Yim{9,(p)ma - lYirni9,ip)ma] (8) 

NFFT in 2-dimensions works on 2-d torus, we have 
thus rewritten equations ([2][8l) into a form (Appendix [B| 
that is suitable to simulate unlcnsed CMB maps at irreg- 
ularly spaced grid points using NFFT. This is possible 
because a band-limited function on a unit sphere can be 
rewritten as a band-limited function on a 2-d torus. In 
order to do this, we have exploited the relation of spin- 
weighted spherical harmonics to Wigncr rotation matri- 
ccs (|A3|) and the factorization of Wigncr rotation matrices 
into two separate rotations (|A5p . 

Using the identities of the spherical triangle, lensed 
temperature anisotropics and polarization in a particu- 
lar direction (9, p) are given by unlensed temperature 
anisotropics and polarization in another direction at the 
last scattering surface. 

f{9,v) ^ T{9',ip') (9) 
Pi9,ip) = exp[-2i(7-a)] P{9\ip') (10) 

The angular coordinates corresponding to the modified 
direction of the photon path {9',p') due to Icnsing arc 
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The Euler angles a, (3 and 7 are defined as. 
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/3 (0 < /? < tt) determines the angle between the direc- 
tions {9,(p) and {9',(p'). a (0 < a < 27r) is the angle 
required to rotate the basis vector 9 = {9, ip) in a right- 
handed sense about h onto the tangent (at n ) to the 
geodesic connecting n and n'; 7 (0 < 7 < 27r) is defined 
in the same manner as a but at n'. 

To compute lensed CMB fields at a particular posi- 
tion on the sphere it is enough to compute the unlensed 
CMB at some other position on the sphere determined 
by the identities of the spherical triangle. The most pop- 
ular pixelization scheme that is used in CMB analysis is 
the HEALPixjs^l pixelization [s^l which is an irregular 
grid on the surface of the unit sphere in {9, (p) coordi- 
nates. Since gravitational lensing remaps the CMB sig- 
nal, the modified angular coordinates due to lensing will 
not, in general, correspond to any other pixel center of 
the HEALPix grid, even if the unlensed CMB is defined 
over HEALPix grid points. Hence, in order to compute 
lensed CMB field on HEALPix grid points, we should be 
able to resample the unlensed CMB at arbitrary posi- 
tions on the sphere. Since remapping on a sphere can be 
recast into remapping on a 2-d torus (Appendix |^ , we 
have used NFFT to compute lensed CMB anisotropics at 
HEALPix grid points. 



IV. SIMULATION OF LENSED CMB MAP 

A. How to simulate a lensed map 

We have seen in the last section that, in the 
Born approximation, gravitational lensing of the CMB 
anisotropics results in a simple resampling of the un- 
lensed anisotropics, with an extra rotation in the case 
of polarization lensing. Let us summarize here the main 
steps of the simulation procedure of lensed CMB maps: 
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• Generate a realization of the (unlensed) CMB har- 
monic coefScients (both temperature and polariza- 
tion) from their (unlensed) power spectra 

• Generate in the same way the harmonic coefficients 
of the lensing potential, or alternatively extract 
them from an N-body simulation 

• Transform the harmonic coefficients of the unlensed 
CMB fields into their 2-d torus Fourier counter- 
parts using equations (jB5l and IB6|) . Also get the 
Fourier coefficients of the displacement field from 
the harmonic coefficients of the lensing potential 
using equation (jBSp . 

• Sample the displacement field at HEALPix centers 
(using eauation lB4l and NFFT), apply this displace- 
ment field to HEALPix pixel centers to get dis- 
placed positions on the sphere (using equations 1111 
fc [T^ . Also compute the extra rotation that will 
be needed for the polarized fields (using equations 
[T3TT4l fc[T5l). 

• Resample the temperature and polarization fields 
at the displaced positions using equations (jBll and 
IB2[) and NFFT, apply the extra rotation to the 
polarized fields. This gives us the simulated Icnsed 
CMB fields, sampled at HEALPix pixel centers. 




Figure 1: (a) A realization of an unlensed CMB temperature 
anisotropics map (nside=1024) that we have obtained using 
NFFT [oversampling factor {a)=2, Convolution length (K) = 
4]. (b) Difference of unlensed CMB maps (nside=1024) that 
we have obtained using NFFT and HEALPix from the same 
harmonic coefficients. 



B. Validation of the method on a known case: 
unlensed maps 

In order to test the part of the algorithm that goes 
from harmonic coefficients of temperature or polarization 
fields to arbitrary real space sampled positions, via 2-d 
torus Fourier modes and NFFT transform, we test the 
method on unlensed temperature or polarization fields, 
sampled at HEALPix centers. Indeed, this is a valid test 
of the method as HEALPix pixel centers are irregularly 
distributed in (9, ip) coordinates. In addition, we can 
directly compare the output of the method to a direct 
resummation of the spherical harmonics decomposition of 
the fields at HEALPix centers by using the fast spherical 
harmonics transforms of the HEALPix package, which 
will serve as a reference. 

In Figure [H we show an (unlensed) realization of 
the CMB temperature anisotropies obtained using our 
method, as well as a map of the difference between our 
method and the HEALPix reference map. Note the dif- 
ference in the color scales. In order to quantify more 
precisely the accuracy of our method, we have computed 
two kinds of error statistics: 



X 



E. 



\X 



NFFT(Oj,(Pj 



HEALPix 



maxn 



X 



HEALPix 



Npix 



X 



NFFT\ 



X 



HEALPi 



ENpix 
.7 = 1 



X 



HEALPi 



where X stands for T, Q, U, dg and d^. 

E^ is the maximum (relative) error for field X, while 
E2 is the relative root mean square error. Table U gives 
the value of these statistics for unlensed CMB tempera- 
ture only. Values of these error norms for the displace- 
ment field and the unlensed CMB polarization fields are 
of the same order of magnitude. 



Table 1: Variation of typical order of magnitude of error norms 
with the convolution length (K) for an unlensed CMB map 
simulated using NFFT. E2 and Eoo are the quadratic norm 
and infinite norm respectively 



Oversampling 


Convolution 


nside 


Imax 


Maximum 


R. M. S 


factor 


length 






error 


error 


{^) 


(K) 






(£00) 


{E2) 


2 


4 


1024 


2048 


~ 10"* 


~ 10"* 


2 


6 


1024 


2048 


~ 10-^^ 


~ 10-^2 


2 


8 


1024 


2048 


~ 10-" 


~ 10-13 



To achieve this accuracy, we have used the Kaiser- 
Bessel window [s^ as the NFFT interpolating func- 
tion. Since the full precomputation of the window func- 
tion at each node in spatial and frequency domains 
requires lots of memory space, we have used a ten- 
sor product form for the multivariate window function, 
that requires only unidimensional precomputations. This 
method uses less memory at the price of some extra mul- 
tiplications [H, [ll] . The accuracy [H, [s^ of our Simula- 



tion can be improved by increasing both the oversamphng 
factor and the convolution length, at the price of extra 
memory consumption and CPU time (Tabic HJ . 



C. Simulation of lensed maps 

We applied our simulation algorithm of lensed CMB 
maps (both temperature and polarization), as described 
in Section IIV Al to 1000 independent realizations with 
HEALPix resolution nside — 1024, and a maximum mul- 
tipole Imax = 2048. In Figure [21 we show one such real- 
ization of a lensed CMB temperature field, as well as the 
difference between the lensed and unlensed fields. 







Figure 2: (Top) A realization of lensed CMB map 
(nside=1024). (Middle) A realization of the amplitude of the 
deflection field (nside=1024). (Bottom) Difference of lensed 
and unlensed CMB maps (nside=1024). These maps are ob- 
tained using NFFT for the oversampling factor {<j)~2 and 
Convolution length (K) = 4. 



Since weak lensing of CMB is a tiny effect at small 
angular scales, we have shown a realization of a small 
portion of the unlensed CMB temperature anisotropics, 
lensed CMB temperature anisotropics, amplitude of de- 
flection field and, the difference of lensed and unlensed 





Figure 3: (Top left) A small portion of a simulated unlensed 
CMB temperature anisotropy map. (Top right). A small por- 
tion of the corresponding lensed CMB temperature anisotropy 
map. (Bottom left) A small portion of the amplitude of the 
simulated deflection field map. (Bottom right) A small por- 
tion of the difference of simulated lensed and unlensed CMB 
maps. These maps are obtained using NFFT for the oversam- 
pling factor {a)— 2 and Convolution length (K) = 4. 



CMB temperature anisotropics in Figure [3] to better il- 
lustrate the lensing effect. Although unlensed and lensed 
CMB temperature anisotropics are indistinguishable to 
the naked eye, the correlation between the deflection field 
and the difference of lensed and unlensed CMB temper- 
ature anisotropics is clearly visible. 

Table [ll] shows the typical CPU time and memory 
that are required to simulate a single realization of un- 
lensed and lensed CMB temperature and polarization, 
with different resolutions. Storage of the window func- 
tion at the grid points both in spatial and frequency do- 
main before computing the Fourier transform consumes 
a fair amount of memory, which ultimately increases the 
overall memory requirement for the simulation of lensed 
CMB maps Hi Hi]. 

Table shows the same, but with different convolu- 
tion lengths. Increase of the convolution length not only 
increases the computational cost of the interpolation part 
of NFFT, but also increases the cost of the precomputa- 
tion of window function and memory requirement as one 
has to compute and store the window function at a larger 
number of grid points in the spatial domain before doing 

nfftP,[3|. 

Figured] shows on the same plots the theoretical power 
spectra C^^, where XY stands for TT, EE, TE, BB re- 



6 



Table II: Variation of CPU time and memory require- 
ments with resolution to simulate CMB maps (unlensed and 
lensed) [oversampling factor (a)— 2, Convolution length (K) — 
4] using NFFT. 



nside 


^max 


CPU 


Memory 






time 


requirement 


256 


512 


1 min 12 sec 


491 MB 


512 


1024 


6 min 8 sec 


1.9 GB 


1024 


2048 


32 min 


7.6 GB 



Table III: Variation of CPU time and memory requirements 
with the convolution length (K) for simulating a realization 
CMB map(unlensed and lensed) [nside= 1024, Imax = 2048] 
using NFFT. 



Oversampling 


Convolution 


CPU 


Memory 


factor 


length 


time 


requirement 




(K) 






2 


4 


32 min 


7.6 GB 


2 


6 


45 min 


8.4 GB 


2 


8 


60 min 


9.1 GB 



spectively, for the lensed and unlensed cases, as predicted 
by CAMB. In the cosmological model we have chosen 
there are no primordial tensors, hence C^^ is entirely 
due to lensing. 

An accurate recovery of this power spectrum from 
lensed polarization maps is therefore a powerful test of 
our simulation method. In Figure [S] we show, on top of 
the lensed theoretical spectra (solid lines), the average 
empirical power spectra computed from the 1000 simula- 
tions (dots). We can see that the agreement is excellent, 
which is remarkable for C^^ as explained above. We 
have ignored the lensed angular power spectrum beyond 
the multipole / = 1700 in the comparison of average em- 
pirical power spectra and theoretical power spectra be- 
cause the accurate computation of the average empiri- 
cal power spectra for the multipolcs / > 1700 requires 
lensed CMB maps simulated from the power spectra of 
unlensed CMB and lensing potential beyond the multi- 
pole I = 2048. 

To have a more quantitative view of the accuracy of 
the method, we show in Figure [H] the relative difference 
between the average empirical power spectra computed 
on the 1000 simulations and the theoretical spectra from 
CAMB, both for the unlensed (red) and lensed (green) 
cases. On each plot, we also show the theoretical r.m.s. 
deviation of the averaged empirical spectra, computed 
neglecting the small lensing-induced non-Gaussianity in 
the lensed cases. Note that this corresponds to a very 
small underestimation of the scatter ([5l|, [s^). Taking 
into account the fact that the averaged power spectra 
are nearly Gaussian distributed (due to the central limit 
theorem), we can assess the presence of possible biases 
in the recovered spectra by computing the reduced 



statistics. 

{Imax - 1) \(C^"^V -I- C^-^C^^l 

1=2 1^1,^;, th ) ^ ^l,t\i '~^l,th 

Here N^iz is the number of independent realizations of 
angular power spectra under consideration. 

Table IV: Reduced statistics for the recovered unlensed 
angular power spectrum 



Angular 


Value of 




power 




P{oo > ZjcY > zxy) 


spectrum 


statistics 




(C'l^aimul) 






,simul 


0.9574 


92% 


.simul 


0.9879 


65% 


^Z ,simul 


0.9901 


62% 



Table V: Reduced statistics for the recovered lensed angu- 
lar power spectrum 



Angular 


Value of 




power 




P(oo > ZjcY > ^xy) 


spectrum 


statistics 




(C'z^imul) 


(4y) 




^Z ,simul 


1.0030 


46% 


f^EE 
^Z ,simul 


0.9928 


58% 


f-^BB 
^Z ,simul 


0.9928 


58% 


pTE 
^Z .simul 


0.9933 


57% 



Table [TVl fc IVl show that the probability of the reduced 
statistics {Z\y) having values greater than the es- 
timated values (zxy) ^^'^ quite large both for unlensed 
and lensed power spectrum. This strengthens our claim 
about the unbiasedness in the simulation of unlensed and 
lensed CMB maps using NFFT. 

V. SUMMARY 

Accurate predictions for the expected CMB 
anisotropics are required for analyzing future CMB 
data sets, which ultimately require accurately simulated 
lensed maps. The most popular pixelization used to 
analyze full-sky CMB maps is the HEALPix pixeliza- 
tion. In order to simulate lensed CMB anisotropics at 
HEALPix grid points wc have to compute unlensed CMB 
anisotropics at irregularly spaced grid points over the 
sphere, determined by the deflection fleld and remapping 
equations. Since remapping on a sphere can be recast 
into remapping on a 2-d torus, we have used the NFFT 
library to compute lensed CMB anisotropics at HEALPix 
grid points and experimented with different settings of 
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Figure 5: Green solid line is the theoretical angular power spectrum Ci_th of lensed CMB, Magenta dots are the average angular 
power spectrum Ci,simui recovered from 1000 realizations of lensed CMB maps (nside=1024 and Imax ~ 2048) 
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Figure 6: Fractional difference of average angular power spectrua recovered from 1000 realizations of CMB maps (nside=1024 
and Imax = 2048) and their corresponding theoretical angular power spectra. Red lines are for unlensed maps and green lines 
are for lensed maps. Black solid lines show the theoretical cosmic variance. 



the accuracy parameters. We have obtained that for a 
nside=1024 map a 10~^ accuracy is easily reached when 
setting the (cr, K) parameters to (2, 4). With our current 
implementation of the method, correspond to a 32 min 
computation on a classical PC configuration. This can 
probably be improved by parallelizing the algorithm. 
Furthermore, the average angular power spectra Cz^simui 
recovered from 1000 realizations of lensed and unlensed 
CMB maps are also found to be consistent with the 
corresponding theoretical ones, Ci^th- This validates our 
simulation of lensed CMB maps. 



Such simulations will be a useful tool for the analysis 
and interpretation of upcoming CMB experiments such 
as PLANCK and ACT. However, they arc not the only 
possible use of this technique. Indeed, the simulation 
of the lensing deflection field can be improved by go- 
ing from the simple Born approximation to ray-tracing 
through dark matter N-Body simulations. Ray-tracing 
faces a similar problem as the simulation of the lens ef- 
fect on CMB maps, i.e. accurately resampling a vector 
field on the sphere. Current state-of-the-art ray-tracing 
algorithms, like (ssj could be made more accurate by us- 
ing the technique described here. 



Appendix A: SPIN s FUNCTIONS ON A SPHERE 
AND 2-D TORUS 

Spin s square-integrable functions s/(^, <p) on a unit 
sphere are conveniently expanded in spin-weighted spher- 
ical harmonics sYim{Oi'f) of same spin [ssl. [sol. [40| . 

sf{e,^) = £ ^ sflm sYl^{0,^) (Al) 
1=0 ni=-l 

with the inverse transform. 

ji,n = f dn j{e,^),YiUe,^) (A2) 

Jn 

These harmonics, with Z G N, m £ Z and max 
(|"T.|,|s|) < form an orthonormal basis for the decompo- 
sition of spin s square-integrable functions on the sphere. 
They are explicitly given in a factorized form in terms of 
the Wigner rotation matrices D^^^,{{p, 9, p), 

sYi^{9,^) = (-iy^^d±lD*J^_^^{^,e,o) (A3) 

With our conventions for the Euler angles [i^, , we 
have, 

Dl „, (^, 9, p) = e-' ™ ^ dl {9) e-* ™' " (A4) 
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These rotation matrices (jA4p basically characterize the 
rotation of spin- weighted spherical harmonics. Decompo- 
sition shown in equation (jA4p is exploited by factoring 
the rotation matrices in two separate rotation matrices 
as follows H, \M 



X D'^„^, (0,|,p+|)(A5) 



Expressing the the Wigner rotation matrices (jA4[) in 
the above manner (|A5[) . equation (|Aip can be rewritten 
as, 



sf{e,^)^ E E 'e'^"''^'"''^ (A6) 

where 



s/77 



E 

Z— ma3:(|m|,|m'|,|s|) 



21 + 1 



sfh 



X d 



I' m ( 2 ) 



exp 



_^(„^ + s)_ (A7) 



The advantage of factoring the rotation matrices in 
this manner is that now the Eulcr angles only occur in 
complex exponentials and we need to evaluate dmm'{0) 
at e = f only [H, [H [H, [H, p . 

Computation of ip) using equation (|A6p may not 
be the most efficient way, but the presence of expo- 
nentials may be exploited such that techniques of fast 
Fourier transform either on irregular or regular grid may 
be used for rapid computation of double summations 
simultaneously. In both cases, the domain of spin s 
function sf{0, ^) must be extended from the sphere, 
{9, if) e [0, tt] X [0, 27r] to the 2-dimensional torus, [9, (p) £ 
[0, 27r] X [0, 27r] using the symmetry syim{2n — 9,Tr + 
if) = (—1)" sYimi9, ip) of spin- weighted spherical harmon- 
ics so that equation (|A6p becomes a complex to complex 
Fourier transform over 2-dimensional torus. The compu- 
tation of sfmm' for (I'Tij, \m'\) < Imax, iuvolvcs perform- 
ing a 1-dimensional summation over a 2-dimensional grid, 
hence it is of order 0{l^^^). 



m^)= E E 

m— — lrriax "f^' — — ^max 



gi{m ip-\-7n' 



(B3) 



{de + id^) {9, If) 

= E E G^™' e'('"'^+™''') (B4) 

and the corresponding Fourier modes are given by, 



E 

l—'max{\rn\ , |m' | ) 



2/ + 1 

47r 



Tlr. 



^ "m' m \^ 2 / ° V 2 / 







) exp 


—im— 

V 2J 



E 

l—7rLax{\7n\,\rn' | ,2) 



2/ + 1 
47r 



{Elm + i Bhn) 



X d 



™'™ V2/ 



dL, 



rn' (-2) \^ 2 



TT exp 



—im— 
2. 



77)/ 



2/ + 1 



y 

^ V 47r 

/— 77iaa:(|m|,|777'|) 
^ "m' 771 







TT" 




exp 


—im— 
V 2J 



Gfnrn' — ^ ^ (~^) 

/^maa:(|??i|,|m'|,l) 



U; + l)(2; + l) 



X dL 



(I 







TT" 


) (2) 


exp 


—im— 
I 2\ 



(B5) 



(B6) 



(B7) 



(B8) 



Appendix B: CMB FIELDS ON 2-D TORUS 

Factoring the rotation matrices in two separate rota- 
tion matrices (|A4[) and extending the domain of CMB, 
lensing potential and defection fields from sphere to 2-d 
torus, equations (|2|3|7|8p can be rewritten as, 

Ti9,ip)= E Tram' e'^^^+^''^ (Bl) 

P{9,ip)^ J2 E e'<"^+"''^ (B2) 
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